The Mitchell River is a river located in Far North Queensland, Australia. The river rises on the Atherton Tableland about 50 kilometres (31 mi) northwest of Cairns, and flows about 750 kilometres (470 mi) northwest across Cape York Peninsula from Mareeba to the Gulf of Carpentaria. we will use this case study to present some functionalities of the prioriactions package.

We started loading libraries:

library(prioritizr) #To use the simulate features distribution
library(prioriactions)
library(raster) #To plot of shapefiles
library(tmap) #To create cool maps
library(scales) #To standardize the value of amount
library(reshape2) #To use the melt function
library(sp) #To use the spplot function
library(viridis) #To use viridis pallete

Preparing data inputs

We divided the whole catchment (71,630 km2 into 2316 sites (i.e., sub-catchments), each one included the portion of river length between two consecutive river connections. We considered four major threats to freshwater fish species in the catchment: water buffalo (Bubalis bubalis), cane toad (Bufo marinus), river flow alteration (caused by impoundments, channels for water extractions and levee banks) and grazing land use. Also, we used the modelled spatial distribution of 45 fish species in the Mitchell river catchment as our conservation features. The distribution of features and threats () will be simulated while the other input files a prioriactions will be loaded directly from the package as follows:

path <- system.file("extdata/input_big/", package = "prioriactions")

pu_data <- data.table::fread(paste0(path,"unitCost_Big.csv"), data.table = FALSE)
features_data <- data.table::fread(paste0(path,"target_Big.csv"), data.table = FALSE)
bound_data <- data.table::fread(paste0(path,"boundary_Big.csv"), data.table = FALSE)
sensitivity_data <- data.table::fread(paste0(path,"sensibility_Big.csv"), data.table = FALSE)

We load the shapefile of the case study also included as part of the package installation:

#reading shapefile
shp_mitchell = raster::shapefile(paste0(path,"Mitchell.shp"))

sp::spplot(shp_mitchell, zcol = "Shape_Area", names.attr = "Area",
      main = "Planning unit areas", col.regions = viridis::viridis(20))

To model the distribution of features and threats, we use the function simulate_species from prioritizr package. Due this simulation is performed on a raster object, we rasterize the shapefile (through the rasterize function of the raster library) and then vectorize it again with the extract function of the same library.

#We set the seed to can replicate the same simulations
set.seed(1)

#We transform the shapefile to raster with a convenient resolution for simulation
ext <- raster::extent(shp_mitchell)
r <- raster::raster(ext, res = 1000) 
r <- raster::rasterize(shp_mitchell, r, field=1)

#Simulation of 45 features distribution
dist_features <- prioritizr::simulate_species(r, 45, model =    RandomFields::RMgauss(scale = 40000))
#> .............................................

#We project the simulated rasters to shapefile
rij_data <- raster::extract(dist_features, shp_mitchell, fun = mean)
colnames(rij_data) <- features_data$id

In this way we obtain spatially auto correlated features distributions. Thus, for example, the distribution of simulated feature 1 is:

shp_mitchell$specie_distribution <- rij_data[,1]

#Plot distribution with tmap library
tmap::tm_shape(shp_mitchell) +
  tmap::tm_fill("specie_distribution", pal = c("white", "seagreen")) +
  tmap::tm_borders(col="black", lwd = 0.5)

Instead, the distribution of simulated threat 1 is:

#Simulation of 4 threats distribution
dist_threats <- prioritizr::simulate_species(r, 4, model =  RandomFields::RMgauss(scale = 40000))
#> ....

#We project the simulated rasters to shapefile
threats_data <- raster::extract(dist_threats, shp_mitchell, fun = mean)
colnames(threats_data) <- 1:4

shp_mitchell$threat_distribution <- threats_data[,1]

#Plot distribution with tmap library
tmap::tm_shape(shp_mitchell) +
  tmap::tm_fill("threat_distribution", pal = c("white", "red4")) +
  tmap::tm_borders(col="black", lwd = 0.5)

1) Model with presence / absence of threats and features

For the first prioritization exercise, we transform the amount data to values 0 and 1 for both features distribution and threats. In the case of features, a threshold of 0.5 was set, that is, if the value of the amount is equal to or greater than 0.5, it is considered that the features is present at the site. Instead, a threat is considered present if it has any non-zero amount value.

#Modifying continuous to binary amount data of features
rij_data_binary <- ifelse(rij_data > 0.5, 1, 0)
shp_mitchell$feature_distribution <- rij_data_binary[,1]

#Modifying continuous to binary amount data of threats
threats_data_binary <- ifelse(threats_data > 0.0, 1, 0)

#Plot example of feature distribution with tmap library
tmap::tm_shape(shp_mitchell) +
  tmap::tm_fill("feature_distribution", pal = c("white", "seagreen"), labels = c("0", "1"), breaks = c(0,1,2)) +
  tmap::tm_borders(col="black", lwd = 0.5)

Then, we create the features and threat distribution inputs in the corresponding format that will be used by the prioriactions package. Also, we modify the features_data to set the target of each feature at 15%.

#Feature distribution data-------------------------------------
#We get the extended matrix of distribution data for both cases, set the specific col names and we are left with only positive values
rij_data_binary_large <- reshape2::melt(rij_data_binary)
colnames(rij_data_binary_large) <- c("pu", "species", "amount")
rij_data_binary_large <- rij_data_binary_large[rij_data_binary_large$amount > 0, ]

#threat distribution data-------------------------------------
threats_data_binary_large <- reshape2::melt(threats_data_binary)
colnames(threats_data_binary_large) <- c("pu", "threats", "amount")
threats_data_binary_large <- threats_data_binary_large[threats_data_binary_large$amount > 0, ]
threats_data_binary_large$cost <- 1
threats_data_binary_large$status <- 0

#features data-------------------------------------
features_data$target <- colSums(rij_data_binary)*0.15

Base model

The base model considers the prioritization of conservation actions using the presence / absence of features and threats. In addition, a target has been set for all features of 15% of their representatives in the basin. The model, for its part, does not consider spatial requirements to the actions (blm and blm_actions parameters equal to 0. Furthermore, there is no attempt to optimize that planning actions occur within the same planning units (i.e., the curve parameter is equal to 1). The solver was configured to stop when a gap of at least 5% is achieved (0% meaning that at least one of the optimal solutions has been found).

input_data.base <- prioriactions::problem(
  pu = pu_data, features = features_data, rij = rij_data_binary_large,
  threats = threats_data_binary_large, sensitivity = sensitivity_data,
  bound = bound_data
)
#> Registered S3 methods overwritten by 'prioriactions':
#>   method                    from      
#>   print.ConservationProblem prioritizr
#>   print.OptimizationProblem prioritizr
#> Correctly loaded inputs
model.base <- prioriactions::min_costs(input_data.base, blm = 0, blm_actions = 0, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.base <- prioriactions::solve(model.base, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 2361 rows, 11580 columns and 106792 nonzeros
#> Model fingerprint: 0x3d759f11
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#> Coefficient statistics:
#>   Matrix range     [3e-01, 4e+00]
#>   Objective range  [1e+00, 1e+00]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+02, 2e+02]
#> Found heuristic solution: objective 3304.0000000
#> Found heuristic solution: objective 2145.0000000
#> Presolve removed 0 rows and 35 columns
#> Presolve time: 0.07s
#> Presolved: 2361 rows, 11545 columns, 106757 nonzeros
#> Variable types: 0 continuous, 11545 integer (11545 binary)
#> 
#> Root relaxation: objective 8.769123e+02, 3052 iterations, 0.19 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0  876.91232    0  694 2145.00000  876.91232  59.1%     -    0s
#> H    0     0                    1397.0000000  876.91232  37.2%     -    0s
#> H    0     0                    1396.0000000  876.91232  37.2%     -    0s
#>      0     0 1029.85774    0  840 1396.00000 1029.85774  26.2%     -    0s
#> H    0     0                    1362.0000000 1029.85774  24.4%     -    0s
#>      0     0 1036.03426    0  885 1362.00000 1036.03426  23.9%     -    1s
#>      0     0 1134.80250    0  739 1362.00000 1134.80250  16.7%     -    1s
#> 
#> Cutting planes:
#>   Cover: 1352
#>   Clique: 10
#>   MIR: 36
#> 
#> Explored 1 nodes (7375 simplex iterations) in 1.38 seconds
#> Thread count was 2 (of 4 available processors)
#> 
#> Solution count 5: 1362 1396 1397 ... 3304
#> 
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.362000000000e+03, best bound 1.135000000000e+03, gap 16.6667%

Note that we have achieved a gap of 0.24% and a objective value of 1265. To obtain the distribution of conservation actions we use the getSolutionActions() function, while that to obtain the vector of sites that are incorporated into the conservation plan we use the getSolutionUnits() function. Note that planning units may be selected just for adding connectivity to the plan without taking action within these sites.

#Getting unit distribution selected
solution_units.base <- solution.base$getSolutionUnits()

#Assign solution to shapefile field to plot it
shp_mitchell$sol.base <-  solution_units.base$solution

#Plot of action 1 distribution
tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("sol.base", pal = c("white", "gray40"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

Whereas, if we want to obtain the actions to be carried out within the sites:

#Getting action distribution
solution_actions.base <- solution.base$getSolutionActions()

#Assign solution to shapefile field to plot it
shp_mitchell$action_1.base <-  solution_actions.base$`1`
shp_mitchell$action_2.base <-  solution_actions.base$`2`
shp_mitchell$action_3.base <-  solution_actions.base$`3`
shp_mitchell$action_4.base <-  solution_actions.base$`4`

#Actions plots
plot_action1.base <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_1.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

plot_action2.base <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_2.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

plot_action3.base <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_3.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

plot_action4.base <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_4.base", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

tmap::tmap_arrange(plot_action1.base , plot_action2.base , plot_action3.base , plot_action4.base)

Or the distribution of the sum of the actions (higher density of actions):

shp_mitchell$sum_actions.base <-  solution_units.base$solution + solution_actions.base$`1` + solution_actions.base$`2` + solution_actions.base$`3` + solution_actions.base$`4`

#Plot som of actions with tmap library
plot.base <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("sum_actions.base", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

plot.base

The results of this base model will be compared with the subsequent models generated.

Model with different curve param

This model differs from the previous one in that it tries to group conservation actions within the selected sites as part of the management plan (through a non-linear relationship in the calculation of benefits). The latter is done by adding a value other than 1 to the curve parameter. Where: (1) indicates that there is a linear relationship between the quotient between the actions carried out with respect to the possible actions to be carried out with respect to the benefit obtained by a characteristic. (2) indicates a quadratic relationship between these values, and (3) a cubic relationship. In simple terms, the cubic relationship further penalizes not performing all conservation actions at a site for a particular feature.

input_data.curve <- prioriactions::problem(
  pu = pu_data, features = features_data, rij = rij_data_binary_large,
  threats = threats_data_binary_large, sensitivity = sensitivity_data,
  bound = bound_data
)
#> Correctly loaded inputs
model.curve <- prioriactions::min_costs(input_data.curve, blm = 0, blm_actions = 0, curve = 3)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.curve <- prioriactions::solve(model.curve, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 301479 rows, 360551 columns and 1153705 nonzeros
#> Model fingerprint: 0x4d5d908f
#> Variable types: 199412 continuous, 161139 integer (161139 binary)
#> Coefficient statistics:
#>   Matrix range     [4e-02, 4e+00]
#>   Objective range  [1e+00, 1e+00]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 2e+02]
#> Presolve removed 49853 rows and 49888 columns (presolve time = 1s) ...
#> Presolve removed 68004 rows and 68039 columns (presolve time = 2s) ...
#> Presolve removed 68004 rows and 74124 columns (presolve time = 3s) ...
#> Presolve removed 68004 rows and 118792 columns (presolve time = 6s) ...
#> Presolve removed 68004 rows and 152707 columns (presolve time = 7s) ...
#> Presolve removed 140554 rows and 152707 columns (presolve time = 7s) ...
#> Presolve removed 140554 rows and 152707 columns (presolve time = 8s) ...
#> Presolve removed 140554 rows and 152707 columns
#> Presolve time: 8.18s
#> Presolved: 160925 rows, 207844 columns, 741013 nonzeros
#> Variable types: 120829 continuous, 87015 integer (87015 binary)
#> 
#> Deterministic concurrent LP optimizer: primal and dual simplex
#> Showing first log only...
#> 
#> 
#> Root simplex log...
#> 
#> Iteration    Objective       Primal Inf.    Dual Inf.      Time
#>        0    8.7000000e+01   5.665120e+03   8.115843e+10     10s
#>    10566    1.6493948e+03   8.977190e+02   3.622366e+10     11s
#>    15871    1.6847355e+03   6.464839e+02   1.832565e+10     12s
#>    21047    1.7383311e+03   3.691404e+02   1.022477e+10     13s
#>    23073    1.7489576e+03   2.919893e+02   1.960043e+10     14s
#>    24993    1.7558028e+03   2.394743e+02   1.106023e+10     15s
#>    27259    1.7580825e+03   1.751924e+02   4.865399e+09     16s
#>    30577    1.7573666e+03   1.397468e+02   4.231776e+10     17s
#>    33064    1.7615242e+03   1.004403e+02   8.089102e+09     18s
#>    35510    1.7729093e+03   6.790741e+01   4.128746e+08     19s
#>    38664    1.7298006e+03   5.834824e+01   5.082061e+10     20s
#>    40759    1.7414736e+03   4.138714e+01   1.519211e+10     21s
#>    43028    1.7471199e+03   1.674000e+01   1.187171e+10     22s
#>    44854    1.7432058e+03   7.907891e+00   2.539410e+09     23s
#>    46708    1.7129775e+03   1.846839e-03   1.070016e+07     24s
#>    46899    4.0241344e+03   0.000000e+00   1.237047e+04     24s
#>    51428    2.6521439e+03   0.000000e+00   3.729628e+04     25s
#>    54201    2.4894892e+03   0.000000e+00   3.926573e+05     26s
#>    57678    2.4029214e+03   0.000000e+00   1.691656e+05     27s
#>    59896    2.3584194e+03   0.000000e+00   1.116342e+04     28s
#>    62449    2.2881074e+03   0.000000e+00   3.687678e+04     29s
#>    64590    2.2361691e+03   0.000000e+00   2.527561e+04     30s
#>    67292    2.1767205e+03   0.000000e+00   5.255631e+04     31s
#>    69685    2.1264155e+03   0.000000e+00   2.306073e+04     32s
#>    72755    2.0664684e+03   0.000000e+00   2.686298e+04     33s
#>    74663    2.0213556e+03   0.000000e+00   6.719329e+04     34s
#>    77022    2.0004158e+03   0.000000e+00   3.404845e+04     35s
#>    79202    1.9625694e+03   0.000000e+00   1.404740e+05     36s
#>    81281    1.9404075e+03   0.000000e+00   6.751236e+04     37s
#>    83370    1.9075073e+03   0.000000e+00   2.233910e+05     38s
#>    85196    1.8781112e+03   0.000000e+00   3.072513e+04     39s
#>    86778    1.8623475e+03   0.000000e+00   5.848692e+04     40s
#>    88779    1.8478615e+03   0.000000e+00   1.581817e+04     41s
#>    91034    1.8202979e+03   0.000000e+00   4.674761e+04     42s
#>    92316    1.8039199e+03   0.000000e+00   7.399480e+04     43s
#>    93842    1.7893068e+03   0.000000e+00   4.079724e+04     44s
#>    95507    1.7741819e+03   0.000000e+00   9.531235e+04     45s
#>    97406    1.7601923e+03   0.000000e+00   1.973123e+05     46s
#>    98853    1.7457871e+03   0.000000e+00   6.493497e+04     47s
#>   100751    1.7274354e+03   0.000000e+00   6.061209e+04     48s
#>   102767    1.7068078e+03   0.000000e+00   8.597032e+04     49s
#>   104235    1.6937678e+03   0.000000e+00   3.242944e+04     50s
#>   105359    1.6845971e+03   0.000000e+00   1.806713e+04     51s
#>   106915    1.6700077e+03   0.000000e+00   1.898088e+04     52s
#>   109377    1.6536675e+03   0.000000e+00   1.706546e+04     53s
#>   110817    1.6459181e+03   0.000000e+00   3.270242e+04     54s
#>   112683    1.6328625e+03   0.000000e+00   9.579376e+04     55s
#>   113847    1.6228827e+03   0.000000e+00   1.734708e+04     56s
#>   115749    1.6072338e+03   0.000000e+00   8.961769e+04     57s
#>   117172    1.5951116e+03   0.000000e+00   2.660406e+04     58s
#>   118800    1.5815913e+03   0.000000e+00   2.249393e+04     59s
#>   120561    1.5666350e+03   0.000000e+00   2.340476e+04     60s
#>   122398    1.5540525e+03   0.000000e+00   5.955468e+04     61s
#>   123646    1.5478251e+03   0.000000e+00   3.784425e+04     62s
#>   124878    1.5411165e+03   0.000000e+00   7.520535e+04     63s
#>   126638    1.5301915e+03   0.000000e+00   3.504995e+04     64s
#>   128054    1.5241642e+03   0.000000e+00   2.780520e+05     65s
#>   129014    1.5205610e+03   0.000000e+00   5.458844e+04     66s
#> Concurrent spin time: 0.00s
#> 
#> Solved with dual simplex
#> 
#> Root relaxation: objective 1.039718e+03, 30576 iterations, 57.01 seconds
#> Total elapsed time = 70.58s
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 1039.71757    0 11890          - 1039.71757      -     -   73s
#>      0     0 1186.30479    0 8871          - 1186.30479      -     -  139s
#>      0     0 1192.02160    0 8333          - 1192.02160      -     -  174s
#>      0     0 1192.62345    0 8579          - 1192.62345      -     -  179s
#>      0     0 1192.72178    0 8479          - 1192.72178      -     -  182s
#>      0     0 1192.73665    0 8435          - 1192.73665      -     -  184s
#>      0     0 1252.04391    0 6723          - 1252.04391      -     -  257s
#>      0     0 1254.79649    0 6584          - 1254.79649      -     -  266s
#>      0     0 1254.97535    0 6618          - 1254.97535      -     -  270s
#>      0     0 1255.01085    0 6658          - 1255.01085      -     -  275s
#>      0     0 1304.65302    0 5183          - 1304.65302      -     -  414s
#>      0     0 1319.82558    0 5053          - 1319.82558      -     -  449s
#>      0     0 1320.37935    0 5093          - 1320.37935      -     -  465s
#>      0     0 1320.49364    0 5086          - 1320.49364      -     -  476s
#>      0     0 1320.53796    0 5055          - 1320.53796      -     -  491s
#>      0     0 1362.32645    0 3870          - 1362.32645      -     -  588s
#> H    0     0                    1896.0000000 1362.32645  28.1%     -  589s
#> H    0     0                    1860.0000000 1362.32645  26.8%     -  589s
#>      0     0 1379.11885    0 3223 1860.00000 1379.11885  25.9%     -  651s
#>      0     0 1380.24570    0 3094 1860.00000 1380.24570  25.8%     -  674s
#>      0     0 1380.31869    0 3156 1860.00000 1380.31869  25.8%     -  688s
#>      0     0 1405.41632    0 2114 1860.00000 1405.41632  24.4%     -  771s
#> H    0     0                    1624.0000000 1405.41632  13.5%     -  774s
#> 
#> Cutting planes:
#>   Gomory: 191
#>   Cover: 38012
#>   Implied bound: 59
#>   Clique: 199
#>   MIR: 5357
#>   Flow cover: 3408
#>   RLT: 3907
#>   Relax-and-lift: 59
#> 
#> Explored 1 nodes (128733 simplex iterations) in 774.63 seconds
#> Thread count was 2 (of 4 available processors)
#> 
#> Solution count 3: 1624 1860 1896 
#> 
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.624000000000e+03, best bound 1.406000000000e+03, gap 13.4236%

Note that the new model increased its dimensions compared to the previous one (from 11580 variables and 2361 constraints to 207844 variables and 160925 constraints), and therefore, increasing its resolution complexity. In this way, we obtain a value of the objective function of 1497 (greater than the previous one) with a quality of 4.48%.

#Getting unit distribution selected
solution_units.curve <- solution.curve$getSolutionUnits()

#Getting action distribution
solution_actions.curve <- solution.curve$getSolutionActions()

#Assign solution to shapefile field to plot it
shp_mitchell$action_1.curve <-  solution_actions.curve$`1`
shp_mitchell$action_2.curve <-  solution_actions.curve$`2`
shp_mitchell$action_3.curve <-  solution_actions.curve$`3`
shp_mitchell$action_4.curve <-  solution_actions.curve$`4`

shp_mitchell$sum_actions.curve <-  solution_units.curve$solution + solution_actions.curve$`1` + solution_actions.curve$`2` + solution_actions.curve$`3` + solution_actions.curve$`4`

#Plot som of actions with tmap library
plot.curve <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("sum_actions.curve", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.curve)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Thus, we can notice that there are a greater number of sites with a higher density of actions with respect to the previous model (for example, a greater number of units painted yellow).

Model with spatial requirements

To add spatial requirements to the model, there are two key parameters: blm which works in the same way as the parameter of the same name in marxan, which tries to minimize fragmentation between the selected planning units (i.e. regardless if conservation actions are carried out within them). Next we will see what happens when this parameter is relevant in the creation of the model.

input_data.spatial_req<- prioriactions::problem(
  pu = pu_data, features = features_data, rij = rij_data_binary_large,
  threats = threats_data_binary_large, sensitivity = sensitivity_data,
  bound = bound_data
)
#> Correctly loaded inputs
model.spatial_req <- prioriactions::min_costs(input_data.spatial_req, blm = 10, blm_actions = 0, curve = 1)
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.spatial_req <- prioriactions::solve(model.spatial_req, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 43557 rows, 25312 columns and 202916 nonzeros
#> Model fingerprint: 0x3fe989f5
#> Variable types: 0 continuous, 25312 integer (25312 binary)
#> Coefficient statistics:
#>   Matrix range     [3e-01, 4e+00]
#>   Objective range  [1e+00, 1e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 2e+02]
#> Found heuristic solution: objective 3304.0000000
#> Presolve removed 0 rows and 35 columns
#> Presolve time: 0.21s
#> Presolved: 43557 rows, 25277 columns, 202881 nonzeros
#> Variable types: 0 continuous, 25277 integer (25277 binary)
#> 
#> Deterministic concurrent LP optimizer: primal and dual simplex
#> Showing first log only...
#> 
#> 
#> Root simplex log...
#> 
#> Iteration    Objective       Primal Inf.    Dual Inf.      Time
#>    31094    3.0518120e+03   0.000000e+00   1.657631e+05      1s
#>    36117    2.2436091e+03   0.000000e+00   2.109409e+05      2s
#>    40599    2.1094886e+03   0.000000e+00   1.042250e+05      3s
#>    44622    1.9300104e+03   0.000000e+00   4.997226e+04      4s
#>    46823    1.8125023e+03   0.000000e+00   5.160657e+06      5s
#>    48643    1.7211745e+03   0.000000e+00   7.747742e+04      6s
#>    50285    1.6811883e+03   0.000000e+00   1.170770e+07      7s
#>    51800    1.6155709e+03   0.000000e+00   2.597011e+05      8s
#>    53093    1.5878591e+03   0.000000e+00   5.304713e+04      9s
#> Concurrent spin time: 0.60s
#> 
#> Solved with dual simplex
#> 
#> Root relaxation: objective 1.094585e+03, 15227 iterations, 9.03 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 1094.58473    0 12761 3304.00000 1094.58473  66.9%     -   10s
#> H    0     0                    3296.0000000 1094.58473  66.8%     -   11s
#>      0     0 1243.82454    0 18610 3296.00000 1243.82454  62.3%     -   29s
#>      0     0 1490.21388    0 13269 3296.00000 1490.21388  54.8%     -   43s
#>      0     0 1490.21469    0 13271 3296.00000 1490.21469  54.8%     -   56s
#>      0     0 1490.21469    0 13272 3296.00000 1490.21469  54.8%     -   70s
#>      0     0 1490.21469    0 13273 3296.00000 1490.21469  54.8%     -   70s
#>      0     0 1490.21469    0 13269 3296.00000 1490.21469  54.8%     -   80s
#> H    0     0                    2414.0927254 1490.21469  38.3%     -   83s
#> H    0     0                    2267.1903810 1490.21469  34.3%     -   84s
#> H    0     0                    2245.1903810 1490.21469  33.6%     -   84s
#> H    0     2                    2235.1903810 1490.21469  33.3%     -   85s
#>      0     2 1490.21469    0 13267 2235.19038 1490.21469  33.3%     -   85s
#>      1     4 1506.90672    1 13296 2235.19038 1490.24408  33.3%  5636   99s
#>      3     6 1520.70894    2 12947 2235.19038 1507.43883  32.6%  7508  115s
#>      5     8 1596.13436    2 13427 2235.19038 1520.76378  32.0%  7228  127s
#>      7    10 1605.31844    3 12311 2235.19038 1520.78909  32.0%  6821  134s
#>      9    12 1530.71185    3 12332 2235.19038 1530.71185  31.5%  5688  211s
#>     11    14 1551.48967    4 12282 2235.19038 1530.71351  31.5%  5743  222s
#>     13    16 1575.02183    5 12700 2235.19038 1530.71351  31.5%  5761  230s
#>     15    18 1622.37087    5 12084 2235.19038 1530.71351  31.5%  5646  238s
#>     17    20 1628.58446    6 11983 2235.19038 1530.71351  31.5%  5651  249s
#>     19    23 1579.58649    6 12678 2235.19038 1530.71351  31.5%  5598  263s
#>     22    25 1620.57003    8 12483 2235.19038 1530.71351  31.5%  5930  275s
#>     24    27 1595.54701    8 12536 2235.19038 1530.71351  31.5%  5901  286s
#> H   26    29                    2233.1903810 1530.71351  31.5%  5689  301s
#> H   27    29                    2225.1903810 1530.71351  31.2%  5872  301s
#>     28    31 1644.45520    9 11236 2225.19038 1530.71351  31.2%  6237  316s
#> H   29    31                    2193.8795220 1530.71351  30.2%  6021  316s
#>     30    33 1626.77311   10 12249 2193.87952 1530.71351  30.2%  6292  327s
#>     32    35 1655.34224   11 12501 2193.87952 1530.71351  30.2%  6254  337s
#>     34    37 1674.91212   11 10028 2193.87952 1530.71351  30.2%  6394  348s
#>     36    40 1701.99427   12 8894 2193.87952 1530.71351  30.2%  6414  359s
#>     39    44 1735.29735   13 10116 2193.87952 1530.71351  30.2%  6320  375s
#>     43    51 1726.77391   14 7611 2193.87952 1530.71351  30.2%  6189  388s
#>     50    55 1764.48171   18 6790 2193.87952 1530.71351  30.2%  5862  400s
#> H   54    61                    2172.8795220 1530.71351  29.6%  5816  412s
#>     60    69 1779.97493   21 5736 2172.87952 1530.71351  29.6%  5655  425s
#>     68    74 1809.64250   25 3476 2172.87952 1530.71351  29.6%  5321  436s
#>     73    87 1821.36263   26 2316 2172.87952 1530.71351  29.6%  5234  447s
#> H   86    99                    1845.0253024 1530.71351  17.0%  4652  454s
#> H   90    99                    1824.1807319 1530.71351  16.1%  4447  454s
#> H   92    99                    1823.1807319 1530.71351  16.0%  4351  455s
#> 
#> Cutting planes:
#>   Zero half: 1
#> 
#> Explored 98 nodes (452968 simplex iterations) in 455.03 seconds
#> Thread count was 2 (of 4 available processors)
#> 
#> Solution count 10: 1823.18 1824.18 1845.03 ... 2267.19
#> 
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.823180731900e+03, best bound 1.530713509025e+03, gap 16.0416%

In the same way as what happened in the previous model, it is naturally inferred that the spatial requirement implies higher costs in the management plan, as well as a model with greater computational complexity to resolve.

#Getting unit distribution selected
solution_units.spatial_req <- solution.spatial_req$getSolutionUnits()

#Getting action distribution
solution_actions.spatial_req<- solution.spatial_req$getSolutionActions()

#Assign solution to shapefile field to plot it
shp_mitchell$action_1.spatial_req <-  solution_actions.spatial_req$`1`
shp_mitchell$action_2.spatial_req <-  solution_actions.spatial_req$`2`
shp_mitchell$action_3.spatial_req <-  solution_actions.spatial_req$`3`
shp_mitchell$action_4.spatial_req <-  solution_actions.spatial_req$`4`

shp_mitchell$sum_actions.spatial_req <-  solution_units.spatial_req$solution + solution_actions.spatial_req$`1` + solution_actions.spatial_req$`2` + solution_actions.spatial_req$`3` + solution_actions.spatial_req$`4`

#Plot som of actions with tmap library
plot.spatial_req <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("sum_actions.spatial_req", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.spatial_req)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

As seen, minimizing fragmentation produced a more closely related conservation plan compared to the base model.

Another parameter to take into account is the blm_actions that adds the minimization of fragmentation between conservation actions, this means that it tries to spatially bring the management actions closer together and therefore, independently achieve a management plan as well related.

input_data.spatial_req_actions<- prioriactions::problem(
  pu = pu_data, features = features_data, rij = rij_data_binary_large,
  threats = threats_data_binary_large, sensitivity = sensitivity_data,
  bound = bound_data
)
#> Correctly loaded inputs
model.spatial_req_actions <- prioriactions::min_costs(input_data.spatial_req_actions, blm = 0, blm_actions = 10, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
solution.spatial_req_actions <- prioriactions::solve(model.spatial_req_actions, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 167145 rows, 66508 columns and 491288 nonzeros
#> Model fingerprint: 0x8086a2c3
#> Variable types: 0 continuous, 66508 integer (66508 binary)
#> Coefficient statistics:
#>   Matrix range     [3e-01, 4e+00]
#>   Objective range  [1e+00, 1e+02]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+00, 2e+02]
#> Found heuristic solution: objective 11548.000000
#> Presolve time: 0.83s
#> Presolved: 167145 rows, 66508 columns, 491288 nonzeros
#> Variable types: 0 continuous, 66508 integer (66508 binary)
#> 
#> Deterministic concurrent LP optimizer: primal and dual simplex
#> Showing first log only...
#> 
#> 
#> Root simplex log...
#> 
#> Iteration    Objective       Primal Inf.    Dual Inf.      Time
#>        0    4.4417447e+04   0.000000e+00   4.631422e+05      2s
#> Warning: Markowitz tolerance tightened to 0.5
#>    61617    2.4718346e+04   1.030990e-03   2.720499e+07      2s
#>    66858    2.2891668e+04   0.000000e+00   3.852743e+05      3s
#>    69133    2.1911720e+04   0.000000e+00   2.224889e+06      4s
#>    71362    2.1302148e+04   0.000000e+00   4.572488e+05      5s
#>    73592    2.0670783e+04   0.000000e+00   1.453689e+06      6s
#>    75648    2.0151453e+04   0.000000e+00   4.117813e+05      7s
#>    77622    1.9833790e+04   0.000000e+00   4.677604e+05      8s
#>    79253    1.9750604e+04   0.000000e+00   4.578073e+05      9s
#>    80145    1.8872755e+04   0.000000e+00   1.105365e+07     10s
#>    81929    1.7169011e+04   0.000000e+00   1.248981e+07     11s
#>    83267    1.6652827e+04   0.000000e+00   9.086708e+05     12s
#>    84605    1.5942891e+04   0.000000e+00   1.345409e+06     13s
#>    86166    1.5443126e+04   0.000000e+00   5.417270e+06     14s
#>    87504    1.5222457e+04   0.000000e+00   2.737745e+07     15s
#>    88620    1.4570497e+04   0.000000e+00   1.905453e+06     16s
#>    89958    1.4107801e+04   0.000000e+00   1.628268e+06     17s
#>    91296    1.3683932e+04   0.000000e+00   2.214303e+06     18s
#>    92857    1.3290566e+04   0.000000e+00   2.393024e+06     19s
#>    94195    1.3000940e+04   0.000000e+00   4.811312e+06     20s
#>    95533    1.2857727e+04   0.000000e+00   1.515921e+06     21s
#>    97094    1.2696943e+04   0.000000e+00   1.366198e+06     22s
#>    98655    1.2491999e+04   0.000000e+00   2.139967e+06     23s
#>    99547    1.2379188e+04   0.000000e+00   2.647626e+06     24s
#>   100885    1.2050869e+04   0.000000e+00   5.010568e+06     25s
#>   102446    1.1910771e+04   0.000000e+00   1.677002e+06     26s
#>   104007    1.1800496e+04   0.000000e+00   1.166538e+06     27s
#>   105791    1.1582139e+04   0.000000e+00   1.531039e+06     28s
#>   107129    1.1310217e+04   0.000000e+00   4.454977e+05     29s
#>   108467    1.1105586e+04   0.000000e+00   9.934501e+05     30s
#>   110028    1.0944070e+04   0.000000e+00   4.081888e+05     31s
#>   111143    1.0860066e+04   0.000000e+00   8.532527e+07     32s
#>   112704    1.0211889e+04   0.000000e+00   5.370109e+06     33s
#>   114042    9.8247835e+03   0.000000e+00   2.867375e+06     34s
#> Concurrent spin time: 0.00s
#> 
#> Solved with dual simplex
#> 
#> Root relaxation: objective 1.289311e+03, 19774 iterations, 33.49 seconds
#> Total elapsed time = 35.10s
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0 1289.31051    0 17015 11548.0000 1289.31051  88.8%     -   40s
#> H    0     0                    10906.955285 1289.31051  88.2%     -   41s
#> H    0     0                    9770.6621484 1289.31051  86.8%     -   41s
#>      0     0 1322.21108    0 28897 9770.66215 1322.21108  86.5%     -   76s
#>      0     0 1691.79293    0 31021 9770.66215 1691.79293  82.7%     -  169s
#> H    0     0                    5843.0029202 1691.79293  71.0%     -  169s
#> H    0     0                    5842.0029202 1691.79293  71.0%     -  169s
#> H    0     0                    5736.9354972 1691.79293  70.5%     -  169s
#>      0     0 1691.79293    0 36204 5736.93550 1691.79293  70.5%     -  225s
#> H    0     0                    5486.5993358 1691.79293  69.2%     -  225s
#> H    0     0                    5464.3631894 1691.79293  69.0%     -  226s
#>      0     0 1691.82283    0 36205 5464.36319 1691.82283  69.0%     -  226s
#>      0     0 1691.82283    0 36205 5464.36319 1691.82283  69.0%     -  304s
#> H    0     0                    5309.6680442 1691.87933  68.1%     -  338s
#> H    0     0                    4213.3417431 1691.87933  59.8%     -  342s
#> H    0     2                    4213.1863136 1691.87933  59.8%     -  343s
#>      0     2 1691.87933    0 36205 4213.18631 1691.87933  59.8%     -  343s
#>      1     4 1709.96480    1 36294 4213.18631 1691.87933  59.8%  3309  365s
#>      3     6 1730.99919    2 36432 4213.18631 1709.96695  59.4%  3502  393s
#>      5     8 1736.65762    2 26842 4213.18631 1721.29090  59.1%  3103  409s
#>      7    10 1764.96275    3 30424 4213.18631 1731.00579  58.9%  3013  421s
#>      9    12 1758.68503    3 26930 4213.18631 1731.00579  58.9%  2860  435s
#>     11    14 1771.40852    4 23138 4213.18631 1731.00579  58.9%  2831  454s
#>     13    16 1793.83908    5 23135 4213.18631 1731.00579  58.9%  3002  465s
#>     15    18 1832.32558    5 22295 4213.18631 1731.00579  58.9%  2861  476s
#>     17    20 1841.19368    6 21344 4213.18631 1731.00579  58.9%  2845  492s
#>     19    22 1810.76823    6 25890 4213.18631 1731.00579  58.9%  2804  504s
#>     21    24 1846.91216    7 22504 4213.18631 1731.00579  58.9%  2701  515s
#>     23    26 1860.50029    7 23807 4213.18631 1731.00579  58.9%  2653  525s
#>     25    29 1847.60307    8 22424 4213.18631 1731.00579  58.9%  2555  532s
#> H   28    31                    3496.0563978 1731.00579  50.5%  2363  555s
#> H   29    31                    3199.6890502 1731.00579  45.9%  2324  555s
#>     30    35 1868.08674   10 22671 3199.68905 1731.00579  45.9%  2461  570s
#>     34    37 1885.18418   10 20161 3199.68905 1731.00579  45.9%  2350  586s
#>     36    41 1894.64458   11 21364 3199.68905 1731.00579  45.9%  2352  606s
#>     40    45 1894.60195   12 19349 3199.68905 1731.00579  45.9%  2306  627s
#>     44    49 1920.74844   14 20776 3199.68905 1731.00579  45.9%  2285  645s
#>     48    55 1927.96817   15 20650 3199.68905 1731.00579  45.9%  2243  660s
#>     54    59 1971.17510   18 20518 3199.68905 1731.00579  45.9%  2144  670s
#> H   58    63                    3097.7023441 1731.00579  44.1%  2045  686s
#> H   61    63                    2538.3060493 1731.00579  31.8%  2002  686s
#>     62    67 2004.85523   21 20189 2538.30605 1731.00579  31.8%  2024  698s
#>     66    72 2021.82057   22 17772 2538.30605 1731.00579  31.8%  1996  711s
#>     71    75 2033.33996   23 20909 2538.30605 1731.00579  31.8%  1980  768s
#>     74    81 2073.08900   24 23025 2538.30605 1731.00579  31.8%  1959  783s
#>     80    88 2005.55142   25 20002 2538.30605 1731.00579  31.8%  1892  797s
#>     87    92 2014.19431   26 20619 2538.30605 1731.00579  31.8%  1826  822s
#> H   88    92                    2301.6723228 1731.00579  24.8%  1805  822s
#>     91   101 2014.51546   27 19985 2301.67232 1731.00579  24.8%  1783  837s
#>    100   110 2069.48068   30 23060 2301.67232 1731.00579  24.8%  1756  858s
#>    109   119 2064.17616   31 23141 2301.67232 1731.00579  24.8%  1724  875s
#>    118   133 2070.25107   36 17861 2301.67232 1731.00579  24.8%  1668  891s
#>    132   144 2098.22708   41 15863 2301.67232 1731.00579  24.8%  1581  915s
#>    143   163 2122.03727   45 12823 2301.67232 1731.00579  24.8%  1552  942s
#> H  154   163                    2277.9698879 1731.00579  24.0%  1499  942s
#>    162   176 2165.89784   50 12909 2277.96989 1731.00579  24.0%  1471  961s
#>    175   186 2201.05396   54 15225 2277.96989 1731.00579  24.0%  1453  986s
#>    189   188 2247.34693   59 17095 2277.96989 1731.00579  24.0%  1437 1013s
#>    203   191 2268.44703   63 13188 2277.96989 1746.19068  23.3%  1387 1034s
#>    214   196 1756.58558    4 35918 2277.96989 1746.19677  23.3%  1358 1061s
#>    219   201 1786.77530    6 35835 2277.96989 1746.19677  23.3%  1374 1097s
#>    224   205 1812.63049    7 33333 2277.96989 1746.19677  23.3%  1394 1125s
#>    228   209 1835.56037    9 32971 2277.96989 1746.19677  23.3%  1413 1159s
#> H  232   202                    2261.9698879 1746.19677  22.8%  1425 1194s
#>    236   209 1875.26386   12 30344 2261.96989 1746.19677  22.8%  1458 1226s
#>    243   216 1938.96010   13 23585 2261.96989 1746.19677  22.8%  1484 1253s
#>    250   228 1959.08680   15 23496 2261.96989 1746.19677  22.8%  1484 1281s
#> H  262   238                    2258.7204477 1746.19677  22.7%  1467 1322s
#>    273   251 1988.74581   18 27565 2258.72045 1746.19677  22.7%  1473 1352s
#>    286   267 2012.31352   21 29597 2258.72045 1746.19677  22.7%  1465 1386s
#>    302   281 2075.94227   22 12576 2258.72045 1746.19677  22.7%  1460 1421s
#>    316   310 2094.87900   27 15761 2258.72045 1746.19677  22.7%  1472 1455s
#>    345   344 2144.34428   28 20344 2258.72045 1746.19677  22.7%  1428 1499s
#> H  349   334                    2242.1120235 1746.19677  22.1%  1421 1499s
#>    381   333 2132.12647   34 19645 2242.11202 1746.19677  22.1%  1347 1537s
#>    394   339 2149.27546   39 22322 2242.11202 1746.19677  22.1%  1339 1634s
#> H  397   254                    2165.5783120 1746.19677  19.4%  1331 1634s
#> 
#> Cutting planes:
#>   Zero half: 1
#> 
#> Explored 401 nodes (594133 simplex iterations) in 1634.78 seconds
#> Thread count was 2 (of 4 available processors)
#> 
#> Solution count 10: 2165.58 2242.11 2258.72 ... 3496.06
#> 
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 2.165578312040e+03, best bound 1.746196772342e+03, gap 19.3658%

The effect of incorporating this parameter (blm_actions) can be seen better in the action distributions:

#Getting unit distribution selected
solution_units.spatial_req_actions <- solution.spatial_req_actions$getSolutionUnits()

#Getting action distribution
solution_actions.spatial_req_actions<- solution.spatial_req_actions$getSolutionActions()

#Assign solution to shapefile field to plot it
shp_mitchell$action_1.spatial_req_actions <-  solution_actions.spatial_req_actions$`1`
shp_mitchell$action_2.spatial_req_actions <-  solution_actions.spatial_req_actions$`2`
shp_mitchell$action_3.spatial_req_actions <-  solution_actions.spatial_req_actions$`3`
shp_mitchell$action_4.spatial_req_actions <-  solution_actions.spatial_req_actions$`4`

#Actions plots
plot_action1.spatial_req_actions <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_1.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

plot_action2.spatial_req_actions <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_2.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

plot_action3.spatial_req_actions <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_3.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

plot_action4.spatial_req_actions <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("action_4.spatial_req_actions", pal = c("white", "dodgerblue4"), labels = c("0", "1"), breaks = c(0,1,2)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

tmap::tmap_arrange(plot_action1.spatial_req_actions , plot_action2.spatial_req_actions , plot_action3.spatial_req_actions , plot_action4.spatial_req_actions)

2) Model with continuous amount of threats and presence / absence of features

An interesting comparison is how the use of continuous threat intensities influences solutions with respect to the use of binary threat presence / absence data. To do this, we change the column amount for the original output from the simulation of threat distribution, without performing a previous filter, i.e.,

#We get the extended matrix of distribution data for both cases, set the specific col names and we are left with only positive values

#threat distribution data-------------------------------------
threats_data_large <- reshape2::melt(threats_data)
colnames(threats_data_large) <- c("pu", "threats", "amount")
threats_data_large <- threats_data_large[threats_data_large$amount > 0, ]
threats_data_large$cost <- 1
threats_data_large$status <- 0

Thus, the new model is build as follow:

input_data.threat_cont <- prioriactions::problem(
  pu = pu_data, features = features_data, rij = rij_data_binary_large,
  threats = threats_data_large, sensitivity = sensitivity_data,
  bound = bound_data
)
#> Correctly loaded inputs
model.threat_cont<- prioriactions::min_costs(input_data.threat_cont, blm = 0, blm_actions = 0, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.threat_cont <- prioriactions::solve(model.threat_cont, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 2361 rows, 11580 columns and 106792 nonzeros
#> Model fingerprint: 0xf3148484
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#> Coefficient statistics:
#>   Matrix range     [3e-02, 4e+00]
#>   Objective range  [1e+00, 1e+00]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+02, 2e+02]
#> Found heuristic solution: objective 3203.0000000
#> Found heuristic solution: objective 1934.0000000
#> Presolve removed 0 rows and 35 columns
#> Presolve time: 0.08s
#> Presolved: 2361 rows, 11545 columns, 106757 nonzeros
#> Variable types: 0 continuous, 11545 integer (11545 binary)
#> 
#> Root relaxation: objective 7.061076e+02, 4161 iterations, 0.28 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0  706.10764    0  596 1934.00000  706.10764  63.5%     -    0s
#> H    0     0                    1142.0000000  706.10764  38.2%     -    0s
#> H    0     0                    1131.0000000  706.10764  37.6%     -    0s
#>      0     0  832.81049    0  742 1131.00000  832.81049  26.4%     -    0s
#> H    0     0                    1102.0000000  832.81049  24.4%     -    1s
#>      0     0  832.84104    0  742 1102.00000  832.84104  24.4%     -    1s
#>      0     0  932.02978    0  621 1102.00000  932.02978  15.4%     -    1s
#> 
#> Cutting planes:
#>   Cover: 1152
#>   Clique: 4
#>   MIR: 24
#> 
#> Explored 1 nodes (7980 simplex iterations) in 1.58 seconds
#> Thread count was 2 (of 4 available processors)
#> 
#> Solution count 5: 1102 1131 1142 ... 3203
#> 
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.102000000000e+03, best bound 9.330000000000e+02, gap 15.3358%
#Getting unit distribution selected
solution_units.threat_cont <- solution.threat_cont$getSolutionUnits()

#Getting action distribution
solution_actions.threat_cont <- solution.threat_cont$getSolutionActions()

#Assign solution to shapefile field to plot it
shp_mitchell$action_1.threat_cont <-  solution_actions.threat_cont$`1`
shp_mitchell$action_2.threat_cont <-  solution_actions.threat_cont$`2`
shp_mitchell$action_3.threat_cont <-  solution_actions.threat_cont$`3`
shp_mitchell$action_4.threat_cont <-  solution_actions.threat_cont$`4`

shp_mitchell$sum_actions.threat_cont <-  solution_units.threat_cont$solution + solution_actions.threat_cont$`1` + solution_actions.threat_cont$`2` + solution_actions.threat_cont$`3` + solution_actions.threat_cont$`4`

#Plot som of actions with tmap library
plot.threat_cont <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("sum_actions.threat_cont", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.threat_cont)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

3) Model with continuous amount of threats and features

The more general model allows us to work with continuous data on intensity of species and threats. For this, the same procedure is carried out on the data done in the previous step, but now with the amount column of the features. Also, it is necessary to update the target column from the features_data file.

#Feature distribution data-------------------------------------

rij_data_large <- reshape2::melt(rij_data)
colnames(rij_data_large) <- c("pu", "species", "amount")
rij_data_large <- rij_data_large[rij_data_large$amount > 0, ]

#features data-------------------------------------
features_data$target <- colSums(rij_data)*0.15

Thus, the new model is build as follow:

input_data.feature_cont <- prioriactions::problem(
  pu = pu_data, features = features_data, rij = rij_data_large,
  threats = threats_data_large, sensitivity = sensitivity_data,
  bound = bound_data
)
#> Correctly loaded inputs
model.feature_cont<- prioriactions::min_costs(input_data.feature_cont, blm = 0, blm_actions = 0, curve = 1)
#> Warning: The blm argument was set to 0, so the boundary data has no effect
#> Warning: The blm_actions argument was set to 0, so the boundary data has no
#> effect
solution.feature_cont <- prioriactions::solve(model.feature_cont, gap_limit = 0.2, verbose = TRUE, output_file = FALSE)
#> Gurobi Optimizer version 9.1.2 build v9.1.2rc0 (linux64)
#> Thread count: 2 physical cores, 4 logical processors, using up to 2 threads
#> Optimize a model with 2361 rows, 11580 columns and 213072 nonzeros
#> Model fingerprint: 0x73cee9ae
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#> Coefficient statistics:
#>   Matrix range     [2e-03, 4e+00]
#>   Objective range  [1e+00, 1e+00]
#>   Bounds range     [1e+00, 1e+00]
#>   RHS range        [1e+02, 2e+02]
#> Found heuristic solution: objective 3302.0000000
#> Found heuristic solution: objective 2299.0000000
#> Presolve time: 0.20s
#> Presolved: 2361 rows, 11580 columns, 213072 nonzeros
#> Variable types: 0 continuous, 11580 integer (11580 binary)
#> 
#> Root relaxation: objective 8.270739e+02, 2992 iterations, 0.17 seconds
#> 
#>     Nodes    |    Current Node    |     Objective Bounds      |     Work
#>  Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
#> 
#>      0     0  827.07394    0  662 2299.00000  827.07394  64.0%     -    0s
#> H    0     0                    1317.0000000  827.07394  37.2%     -    0s
#>      0     0  959.14938    0  789 1317.00000  959.14938  27.2%     -    1s
#> H    0     0                    1289.0000000  959.14938  25.6%     -    1s
#>      0     0  959.25344    0  790 1289.00000  959.25344  25.6%     -    1s
#>      0     0  959.25585    0  790 1289.00000  959.25585  25.6%     -    1s
#>      0     0 1054.37015    0  788 1289.00000 1054.37015  18.2%     -    2s
#> 
#> Cutting planes:
#>   Cover: 1323
#>   Clique: 2
#>   MIR: 39
#> 
#> Explored 1 nodes (6483 simplex iterations) in 2.29 seconds
#> Thread count was 2 (of 4 available processors)
#> 
#> Solution count 4: 1289 1317 2299 3302 
#> 
#> Optimal solution found (tolerance 2.00e-01)
#> Best objective 1.289000000000e+03, best bound 1.055000000000e+03, gap 18.1536%
#Getting unit distribution selected
solution_units.feature_cont <- solution.feature_cont$getSolutionUnits()

#Getting action distribution
solution_actions.feature_cont <- solution.feature_cont$getSolutionActions()

#Assign solution to shapefile field to plot it
shp_mitchell$action_1.feature_cont <-  solution_actions.feature_cont$`1`
shp_mitchell$action_2.feature_cont <-  solution_actions.feature_cont$`2`
shp_mitchell$action_3.feature_cont <-  solution_actions.feature_cont$`3`
shp_mitchell$action_4.feature_cont <-  solution_actions.feature_cont$`4`

shp_mitchell$sum_actions.feature_cont <-  solution_units.feature_cont$solution + solution_actions.feature_cont$`1` + solution_actions.feature_cont$`2` + solution_actions.feature_cont$`3` + solution_actions.feature_cont$`4`

#Plot sum of actions with tmap library
plot.feature_cont <- tmap::tm_shape(shp_mitchell) + 
  tmap::tm_fill("sum_actions.feature_cont", palette="viridis", labels = c("do nothing", "pu selected without actions", "one action", "two actions", "three actions", "four actions"), breaks = c(0,1,2,3,4,5,6)) + 
  tmap::tm_borders(col="black", lwd = 0.5)

#Comparative with base model solution
tmap::tmap_arrange(plot.base, plot.feature_cont)
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
#> Some legend labels were too wide. These labels have been resized to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.